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Abstract. A coagulation-decoagulation model is introduced on a chain of length L 
with open boundary. The model consists of one species of particles which diffuse, 
coagulate and decoagulate preferentially in the leftward direction. They are also 
injected and extracted from the left boundary with different rates. We will show 
that on a specific plane in the space of parameters, the steady state weights can be 
calculated exactly using a matrix product method. The model exhibits a first-order 
phase transition between a low-density and a high-density phase. The density profile 
of the particles in each phase is obtained both analytically and using the Monte Carlo 
Simulation. The two-point density-density correlation function in each phase has also 
been calculated. By applying the Yang-Lee theory we can predict the same phase 
diagram for the model. This model is further evidence for the applicability of the 
Yang-Lee theory in the non-equilibrium statistical mechanics context. 
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1. Introduction 

It is known that one- dimensional driven diffusive systems can exhibit many interesting 
collective phenomena such as jamming and spontaneous symmetry breaking 
Applications of such models are divers and include the kinetics of biopolymerization 
[Sj and traffic flow |E], but to give just a few examples. 

A powerful framework for studying equilibrium phase transitions in classical statistical 
mechanics is the Yang-Lee theory of phase transition E] • Here we shall briefly explain 
the key points of this theory. Consider a simple model consisting of L spins in thermal 
contact with a heat reservoir. The hallmark of a phase transition in this system is the 
appearance of nonanalyticities in its free energy which can be defined in terms of the 
partition function of the system as / = -^lnZ^. In the thermal equilibrium the 
partition function (which is a function of temperature T in our example) determines the 
statistical properties of the system. However, since the partition function has no real 
and positive zeros in the complex-T plane, there is no scope for a phase transition in a 
finite system. As we increase the number of spins L to infinity, the partition function 
zeros might accumulate towards a point To on the real axis so there is the possibility 
of a phase transition at this point. The density of zeros near the accumulation point 
determines the order of the phase transition. At a first-order phase transition there is 
a nonzero density of zeros at the critical point whereas the density of zeros decays to 
zero at a continuous transition. Although we explained the theory of partition function 
zeros with reference to a system described by a canonical partition function, we should 
note that the theory still holds when one is working in grand canonical ensembles. In 
this case the fugacity (or chemical potential) plays similar role as the temperature T. 
More generally one can look for the zeros of the partition function in the complex plane 
of any of its intensive variables [7J |H] . 

The lack of such a theory in the non-equilibrium statistical mechanics context has 
recently encouraged people to examine whether a generalized version of it can still work 
for out-of-equilibrium systems. Quite similar to the equilibrium statistical mechanics one 
can define a generalized partition function as sum over the steady state weights. These 
weights can be calculated using different methods such as matrix product formalism in 
which the stationary probability distribution function is written in terms of the products 
of non-commuting operators ^3] • Then we look at the zeros of the partition function in 
the thermodynamic limit in the complex plane of any intensive variable of the model. 
The zeros might divide the complex plane into several regions where the free energy of 
the system is analytic. Each region corresponds to a different phase and as we explained 
above, the order of transition depends on how the zeros accumulate towards the positive 
real axis. Several out of equilibrium models have been proposed and studied using the 
generalized version of the Yang-Lee theory |H]-[I2], and it seems that the results obtained 
from application of this theory accord themselves well with the known results obtained 
from other approaches. 

In this paper we will introduce a one-dimensional reaction-diffusion model with open 
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boundary and apply the Yang-Lee theory thereto. This model has already been studied 
with reflecting boundaries ^2JE3 while the open boundary problem has not been solved. 
We shall investigate the effects of the open boundary condition on the phase diagram and 
also the phase transition points of this model. Our work will also be another example 
(with more complicated reaction rules) that shows the applicability of the Yang-Lee 
theory to study the critical behaviors of the out of equilibrium systems. 
In chapter (2) we define the model. The chapter (3) is devoted to the application of the 
matrix product formalism in order to find the steady states weights and therefore the 
partition function of the model. Using this formalism allows us to find the stationary 
weights of the model exactly and calculate many interesting quantities such as the 
density profile of particles on the chain and also the correlation functions. It turns out 
that our model has a first-order phase transition between a low-density and a high- 
density phase provided that there is a constraint on the reaction rates. In chapter (4) 
we shall apply the Yang-Lee theory to our model and calculate the grand canonical 
partition function of the system and its line of zeros in the thermodynamic limit. The 
line of zeros is a circle and the density of zeros is constant over it all. In the Yang-Lee 
theory language this signifies a first-order phase transition. Finally in section (5) we 
shall conclude and explain the generalizations. 

2. The Model 

Consider a one- dimensional lattice of length L. From the left boundary classical particles 
are injected with rate a if the target site is empty. They can also be extracted from 
there with rate /?, provided that it is occupied by a particle. In the bulk of the chain 
particles diffuse to the left and right. When two of them meet, they can merge into 
a single particle. Similarly, a single particle can split into two particles. All of these 
processes take place preferentially to the leftward direction. Specifically, the reaction 
rules in the bulk of the chain are 



diffusion to the left: 


+ 


+ A + 


with rate 
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A + $- 
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coagulation at the right: 
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decoagulation to the right: 




* A + A 


with rate 
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(1) 



in which A and represent an occupied and an empty site respectively. At the left 
boundary we have 

injection at the first site: — > A with rate a , , 

extraction at the first site: A — > with rate (3. 
It is seen that (JIJ consists of transitions which involve only the configurations of nearest 
neighboring sites. For q > 1 the particles have a tendency to move in the leftward 
direction. In the following section we will use the Matrix Product Formalism to find the 
steady states weights of this model. 
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3. Application of the Matrix Product Formalism 



The time evolution operator of our model can be written as 



H 



3=1 



in which 
and 



M L) = / i (L) ®x 0(L - 1) 



where 1 is a 2 x 2 identity matrix, /i is a 4 x 4 matrix for the bulk interactions and 
hS L > is a 2 x 2 matrix for particle input and output at the left boundary. In a basis 
(00, %A, A$, AA) the bulk evolution operator h and read 



/ 



V o 





(A + l)g 
Q 

Aq 











—a 
a 



(3 



-(A + ^g- 1 q 

Aq^ 1 —q + q^ 1 j 

In the following we show that the stationary probability distribution of any configuration 
of our model can be calculated exactly using the matrix product formalism (MPF)[13J. 
Let us first briefly review the MPF. According to this approach, for models with 
Hamiltonian H = YlfZi + h[ L ^ + the stationary probability distribution P(C) 

of any configuration C is assumed to be of the form 

P(C) = ±-{w\ flinD + (l - n )E)\v) (4) 

(Tj = if the site i is empty and = 1 if it is occupied by a particle) with the following 
property 

HP(C) = 0. 

The factor Zl in (jlj) is a normalization factor and can easily be obtained using the 
normalization condition Yle = The matrices D and E are square matrices and 
stand for the presence of a particle and an empty site. These matrices beside the vectors 
(W\ and \V) satisfy the following algebra 



(3) 



h 



E 
D 



E 
D 



E 
D 



E 
D 



E 
D 



E 
D 



< W\h^ 



E 
D 



< W\ 



E 
D 



(5) 
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The auxiliary matrices E and D are also square matrices. For the present model with 
= o, one can easily see from Q and (j3J) that the corresponding algebra is 

[C,C\ = [E,E]=Q 

EC-EC={q + qA + q- l )EC - q(l + A)E 2 - q~ l C 2 (6) 
CE-CE= {q- 1 + q- l A + q)CE - q-\l + A)E 2 - qC 2 . 

with 

(W\((a + {3)E + E- (3C) = (W\C = 

(7) 

E\V) = C\V) = 

in which we have defined C := D + E and C := D + E. Traditionally, one should find 
either a representation for the algebra © and ©/or calculate the stationary weights 
P| directly from the algebra without using any representations. In the following we will 
show that a finite-dimensional representation of the algebra (jOJ and (JZJ) exists under 
special constraint on the parameters a, (3, q and A. One can easily check that for 
q 2 ^ 1 + A the following matrices and vectors 



1+A \ = „ „ ( 1 A \ =, / ^-±A --A 



C= ' , ,C = 0, J E?= I _ 2 ,£ 
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q 1 7 ' ' \ if ' 



(8) 



provide a two-dimensional representation for the algebra (jOJ and ((Zj) provided that 
a = (g^ 1 — q + (3) A. For all rates to be positive we assume q > 1 and (3 > q — q^ 1 . For 
g 2 = 1 + A, C is not diagonalizable and we find another representation 



(9 2 ~1) 2 q 2 -l X 

c i i . r o. £• i -r ;').£=[ «j g 




),C = 0,E= ( 




A N 
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(9) 



In (jHJ) and (jHJ), A is a free parameters. Using the representations (JBJ) and © one can 
easily calculate the density profile of the particles in the stationary state and also their 
density-density correlations. The mean density of particle on the chain at site i, pi, can 
be written in terms of the matrices C and E and the vectors \V) and (W\ 

(W\C*-\C - E)C L -*\V) 

{pl) = W\cW) ' (10) 

By using (jHJ) for q 2 ^ 1 + A and after some algebra we find 

(Pi) = (11) 
A(l + A)~ 1 (l - q 2 + g/?)(A(l + A) L q 2i - (q 2 - 1)(1 + A)y L )g" 2i+1 
qA(q 2 - l)(q 2L - (1 + A) L ) + (3(Aq 2 (l + A) L - q 2L (q 2 -!)(! + A)) ' 
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Two different phases are apparent: q 2 > 1 + A and q 2 < 1 + A. In the thermodynamic 
limit L — > oo we find two different expressions for the mean density of particles at site 
i. In the first phase, which is called the low-density phase, we obtain 

lp.\ — g 2 A(g- 1 -g + ^) /€ , ( ) 

W (l + A)(/3(l + A)-gA) e ° ^ >i + A - ^ 

As can be seen, the density of particles at the left boundary has the largest value and 
will quickly drop to zero in the bulk of the chain with characteristic length £ where 
= | In I- F° r the second phase, the high-density phase, we obtain 
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Figure 1. Monte Carlo simulation results for the density profile of the particles on a 
chain of length L = 100 for q — 1.5, (3 = 0.9 and different values of A. The input rate 
a is given by a = (q^ 1 — q + /3)A. 

where = \ l n t+a I - m high-density phase the density of particles is nearly constant 
in the bulk of the chain equal to and drops to zero near the right boundary. For 
q 2 = 1 + A we use (jHJ) to find the following expression for the mean density of particles 
at site % in the thermodynamic limit 

(p(x)) = TTa (1 " x) for ^ 2 = 1 + A ( 14 ) 

in which x :— |- and < % < 1. In Figure (1) we have plotted the density profile of the 
particles on a chain of length L = 100 for different values of A obtained from the Monte 
carlo simulation. One can readily check that the Monte Carlo results accord with the 
exact analytical results (jl2j) -(|14j ) . In order to see the phase transition more clearly let 
us define the density of particles averaged over the entire system as an order parameter 

P:=^E(P,)- (15) 
i=i 

Plot of (fTK|) as a function of A for q — 3, (3 — 3 and two different values of system 
length L = 200,2000 is shown in the Figure (2). It can be seen that as the length of 
the chain approaches the thermodynamic limit L — ► oo, a discontinuity appears in p 
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Figure 2. Plot of (fH)j) as a function of A for systems of length L = 200 and L = 2000 
with q = 3 and [3 = 3. In the thermodynamic limit L — > oo a first-order phase 
transition occurs at A c = q 2 — 1 = 8. 



at the transition point A c = q 2 — 1 which illustrates a first-order phase transition in 
the system. Using the matrix representation of the algebra © and (j7J) we can also 
calculate two-point density correlation functions of this model in each phase exactly. 
The connected two-point density- density correlation function is defined as 

(piPj)c = (PiPj) ~ (Pi)(Pj) for i < j. (16) 

As we saw in (110 J) , the density of particles on the chain at an arbitrary site i can be 
written in terms of the operators C and E and the boundary vectors (W\ and \ V). The 
expression (pipj) can also be written in terms of the same operators and the vectors as 

. (W\C^(C - E)C^-\C - E)C L ~i\V) 

M = W\cW) • (17) 

By using the matrix representation of the algebra we find 



qA 2 (q 2 -q(3- 1) A(l + A) 



L~2 



{ ^ ^W + d Z-A^-i ^-"(1 + A^)] (18) 

in which Zl = (W\C \V) is given by (|2U|) . Now we look for the thermodynamic limit 
of (fTB^l in each phase. Using (fTTj) . (fTB^) and (fTSj) we obtain the following expressions for 
the thermodynamic limit of the connected two-point correlation function of particles 

(PiPj)c = (YTA~ {Pi)){Pj) - (19) 

One should use (|T2|) or ()13|) for (p^ in each phase. This expression is valid for both 
q 2 > 1 + A and q 2 < 1 + A and has also a simple explanation; namely, in the low-density 
phase the connected two-point correlation function (fTB^) is non-zero only in the close 
vicinity of the left boundary where i < j ~ £. For j 3> £, (piPj)c = 0. Thus, in the 
whole region £ j < L, particles do not feel any correlations. In the high-density phase 
both i and j in (j!6)l should be near the last site of the chain far from the left boundary, 
otherwise (piPj) c = 0. 
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4. Application of the Yang-Lee Theory 

Now we will examine the Yang-Lee theory. We define the partition function of the 
system Z L as the sum of the stationary states weights given by (jH) 



Z L = Y.P(C)= E (W\l[(r l D+(l-r l )E)\V) = (W\C L \V). 



{n=o,i} 



i=l 



Using the matrix representation of the operator C and the vectors \V) and (W\ in 
one obtains 

Aq(q 2 - q/3 - 1)(1 + A) L 



g 2L 4 



(g 2 -i)(/3 + /3A- g A) • (20) 

We look for the zeros of our partition function ()20j) in the complex-g plane at fixed A 
and (3 in the thermodynamic limit L — > oo. In Figure (3) we have plotted the numerical 
estimates of these roots for A = 8, (3 = 3 and L = 100. It can be seen that the roots 
lie on a circle of radius q c = a/1 + A = 3, as we had predicted above. The roots also 
intersect the real-g axis at an angle | which is the sign of first-order phase transition. 
Beside the numerical estimates, one can obtain the line of zeros and also their density 



Im q 



-3-2-10 1 2 3 



Re q 

Figure 3. Plot of the numerical estimates for the roots of l|2()ll in the complex q plane 
for L = 100, A = 8 and (3 = 3. 

near the positive real q axis analytically. Following |H] we define the extensive part 
of the free energy as 



g(q,A,/3) = lim - \nZ L 
in which is given by (|20j). Now the line of zeros can be obtained from 
Re gi = Re g 2 . 



(21) 



(22) 



Here g\ and g 2 are the thermodynamic limits of g on the left and the right hand side 
of the phase transition point. By using (j2*U|) . (|2~T|) and (J2*2*|) and after some calculations 
one finds the line of zeros 



2 i 2 

x + y 



A). 



(23) 
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In (|23|) the parameters x and y are defined as the real and the imaginary parts of q 
respectively. This result is in quite agreement with the numerical estimate of zeros 
given in Figure (3). In order to find the density of zeros /i near the positive real q axis, 
we use [7J E] 

2nti(8,A,P) = -^Im(g 1 -g 2 ) (24) 

where s is the arc length of the line of zeros which is measured from the critical point and 
increases along with the line of zeros. The functions g\ and gi have the same definitions 
as mentioned above. After some algebra we find \i = n ^ +A ■ It is seen that the density 
of zeros is constant everywhere on the circle beside in the vicinity of the positive real q 
axis; therefore, it confirms the existence of a first-order phase transition in our model. 



5. Concluding Remarks 

In this paper we introduced a one-dimensional out-of-equilibrium model consisting of 
one class of particles which diffuse, coagulate and decoagulate on an open chain. The 
left boundary of the chain is assumed to be open so that the particles can enter or leave 
the system. Exact calculations using the MPF show that on a specific plane in the 
parameters space a first-order phase transition takes place between a low-density and a 
high-density phase. The density profile of the particles in each phase is obtained exactly. 
In the low-density phase the density profile of particles far from the left boundary drops 
exponentially to zero with the length scale £ = | In y^aI -1 - We should mention that in 
the stationary state there exists only one characteristic length scale £, while the same 
model without injection and extraction of particles has three different length scales 
|14|llf)|. In the high-density phase the density of particles is nearly constant in the bulk; 
however, it drops to zero near the right boundary with the same length scale £. It is 
shown that the steady state of our model can be written in terms of the superposition of 
shocks [16J. In order to answer the question whether we can apply the classical Yang-Lee 
theory to the out-of-equilibrium systems such as the one introduced here, we have seen 
that the application of this theory gives same results obtained from the MPF approach. 
We calculated the line of partition function zeros exactly. It intersects the positive real 
q axis at right angles at q c = yl + A. The density of zeros is uniform at the cross point, 
which means a first-order phase transition at critical point. Our model, which is more 
complex than those studied before, is another example for the generality of the classical 
Yang-Lee picture of phase transitions in one- dimensional out-of-equilibrium models. It 
can be shown that the MPF also works when the particles are allowed to leave the chain 
from the right boundary; however, there should be still a constraint on the reaction 
rates. The whole phase diagram of this model which can be studied at the mean field 
level or using the Monte-Carlo simulation, is still an open problem. 
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